function [Xfunc_mean_parts] = compute_mean_Xfunc(Xparts, myfunc)
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% This function computes 1/N * myfunc(X_i) 
	% where N is the number of observations (dim1_ii),
	% and myfunc is function passed as argument.
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Inputs:
	% % Xparts:				cell(NumParts,1)
	%   {ii}:					object
	%		.X:						dim1_ii x NumX
	%		.X_FEs:					dim1_ii x NumX_FEs
	% myfunc:				function
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Outputs:
	% Xfunc_mean_parts:		cell(NumParts,1)
	%	{ii}:					object:
	%		.X:						1 x NumX
	%		.X_FEs:					cell(NumX_FEs,1)
	%			{ff}:					NumFE_vals_ff x 1
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
	NumParts = length(Xparts);
	Xfunc_mean_parts = cell(NumParts,1);
	for pp = 1:NumParts
		Xpart_pp = Xparts{pp};
		Xfunc_mean_parts{pp}.X = mean(myfunc(Xpart_pp.X),1); % 1 x NumX
		
		[dim1_pp, NumX_FEs] = size(Xpart_pp.X_FEs);
		Xfunc_mean_parts{pp}.X_FEs = cell(NumX_FEs,1);
		for ff = 1:NumX_FEs
			NumFE_vals_ff = Xpart_pp.NumX_FE_vals(ff);
			% If I expanded FEs into dummies, I would have a number of zeros and a number of ones, out of the dim1_pp observations.
			% Count the relative frequency of ones and the relative frequency of zeros.
			relFreqOnes = accumarray(Xpart_pp.X_FEs(:,ff), ones(dim1_pp, 1), [NumFE_vals_ff 1], @sum)./dim1_pp; % NumFE_vals_ff x 1
			relFreqZeros = 1 - relFreqOnes;  % NumFE_vals_ff x 1
			% Evaluate myfunc(1) and myfunc(0) and multiply by relative frequencies, to get average[myfunc(x)], where x is 0 or 1.
			Xfunc_mean_parts{pp}.X_FEs{ff} = relFreqOnes*myfunc(1) + relFreqZeros*myfunc(0);
		end
	end
end
